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j  Abstract 

\  This  paper  presents  approximations  for  the  rms  error  of  the  maximum 
likelihood  estimator  of  the  direction  of  a  plane  wave  incident  on  a 
random  array.  The  sensor  locations  are  assumed  to  be  realizations  of 
independent,  identically  distributed  random  vectors.  The  second  part  of 
the  paper  presents  an  asymptotically  unbiased  estimator  of  the  noise 


wavenumber  spectrum  from  random  array  data.  J. 
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Introduction 


Sonobuoy  fields  are  used  to  detect  submarines.  Thorn,  Booth,  and 
Lockwood*  have  proposed  that  the  signals  from  randomly  deployed  sonobuoys 
be  coherently  combined  to  make  acoustic  measurements.  They  present  the 
expected  value  and  variance  of  the  pattern  function,  and  the  distribution 
of  the  directivity  index  of  a  three-dimensional  random  array.  In  their 
model,  the  sensor  locations  are  observed  realizations  of  random  variables 
that  may  be  correlated  and  have  different  distributions.  They  define  an 
array  to  be  totally  random  if  the  sensor  locations  are  realizations  of 
independent,  identically  distributed  random  variables.  Several 
stochastic  properties  of  the  sldelobe  pattern  of  a  totally  random  array 
are  given  by  Steinberg. ^ 

The  ratio  of  the  peak  sldelobe  to  the  main  lobe  and  the  directivity 
index  of  an  array  system  are  measures  of  its  ability  to  perform  its 
tasks.  The  generic  signal  processing  tasks  of  an  array  system  are: 

1)  detecting  and  estimating  parameters  of  coherent  wave  signals  that 
impinge  on  the  array;  2)  resolving  multiple  wave  signals;  3)  estimating 
range,  bearing,  or  velocity  of  a  source  that  generates  the  detected 
signal;  and  4)  estimating  the  frequency-wavenumber  spectrum  of  the 
ambient  noise  field.  This  description  of  system  tasks  emphasizes  the 
statistical  nature  of  the  problem  of  measuring  performance,  especially 
for  random  arrays. 
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This  paper  presents  approximations  for  the  mean-square  error  of  the 
maximum  likelihood  estimator  of  the  bearing  of  a  plane  wave  impinging  on 
a  random  array  from  a  distant  source.  The  second  part  deals  with 
estimating  the  ambient  noise's  wavenumber  spectrum. 

1.  Random  Planar  Arrays 

Consider  a  planar  array  of  M  sensors  where  the  sensor  locations 
{(xk,yk)}  are  realizations  of  independent,  identically  distributed  random 
variables  {(Xk,Yk)}.  Assume  for  simplicity  that  the  signal  is  a  single 
frequency  plane  wave  plus  stationary,  zero  mean,  Gaussian  noise.  Let  0O 
denote  the  wave's  direction  of  arrival  with  respect  to  the  x  axis.  This 
angle  is  the  source  bearing  if  the  medium  is  horizontally  homogeneous. 

Let  u>0,  X0,  and  A  denote  the  wave's  frequency,  wavelength,  and  complex 
amplitude,  respectively.  The  signal  at  the  kth  sensor  is 

s(t,xk,yk)  -  AexpIKujQt-K^-iCyyjj)]  +e(t,xk,yk),  (1) 

where  kx  **  (2ir/Ao)cos0o  and  iCy  ■  (2TT/Xo)sin0o  are  the  x  and  y  components 
of  is  the  wavenumber,  and  e(t,xk,yk)  is  a  realization  of  the  noise 
field. 

The  correspondence  between  beamforming  and  frequency-wavenumber 

processing,  and  an  approximation  to  the  maximum  likelihood  (ML)  estimator 

i  M  _  9 

of  0O  have  been  presented  in  a  previous  paper. J  if  pEk«i(xk-x)^  and 
M  _2 

pZk-l'Yk-y)*  are  large,  where  p  is  the  power  signal-to-noise  ratio  in  a 

—  1  M  A 

narrow  band  about  u>0  and  x  ■  M  LevinH  shows  that  the  root 


mean-square  errors  of  the  ML  estimators  of  kx  and  Ky  are  approximately 


-  M  _  "1/2 

rmse  <x  -  [2p  \  (x^-x)2] 
k-1 


M  _  "1/2 

rmse  <y  *  [2p  ^  (yk"y)2) 
k-1 


(2) 


*  M  _  _ 

Moreover,  the  covariance  is  E(kx  -  *x)(iCy  -  Ky)  *  (2pElcm^(xjt-x) (y^-y)  )" 

These  expected  values  are  conditional  on  a  realized  array  geometry,  i.e 

they  are  ex-post  the  deployment  of  the  array. 

To  approximate  these  errors,  assume  that  M  Is  large.  Since  the 

sensors  must  lie  in  some  closed  and  bounded  set,  the  random  variables 

(Xk.Yfc)  are  bounded.  Thus  the  central  limit  theorem  Implies  that 

M  M 

M~lj:k-l(xk“*>2  -  ox2  +  0p(M“1/2)  and  -  oy2  +  0p(M-i'2)f 

where  ox2  and  Oy2  are  the  variances  of  Xk  and  respectively,  and 

Op(M~l/2)  means  that  for  any  e>0,  there  Is  a  Be>0  such  that  the  error  1 

bounded  by  B€M"l/2  with  probability  1-e .  Thus  the  rms  errors 

A  A 

of  kx  and  Ky  are  approximately 


rmse  kx  *  (2pM)"l/20“l 

x 

(3) 

A 

rmse  Ky  *  (2pM)_1 /2a"l 

y 


for  large  M.  The  estimators  are  approximately  uncorrelated  If  the 
coordinate  system  is  rotated  to  make  the  covariance  oTO  »  0  after 
rotation. 


it 
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The  maximum  likelihood  estimator  of  the  bearing  Is 

A  A  A  A  A 

0O  ■  tan~^(<y/<x)  radians.  The  linear  approximation  of  tan~^(Ky/tcx)  - 
tan”l(Ky/Kx)  Is 


-2  -1  * 
(l+Ky^Xv  rMiCv 


-2 


y  y' 


y*x 


(kx-*x)]. 


(4) 


Since  k  and  tc  are  approximately  uncorrelated  If  =  0,  It  follows 

a  y  Ay 

from  (3)  and  (4)  that  when  pMox2  and  pMoy2  are  large, 

A 

E(0O-6O)2  »  (Xq/2it )2(2pM)“^(o”^sin^0o  +  o~^cos^0o).  (5) 

x  y 

Thus  If  ox  ■  Oy  =  a,  then  from  (5) 

A  1/2 

rmse0o  3  XQ(2pM)“  (2ira)_1rads.  (6) 

For  example,  let  a/XQ  «■  12,  M  -  90,  and  p  »  l/4(-6  dB) .  Then  from  (6), 
rmse0o  ■  0.11°  (1.98  x  10“-*  rads).  If  o/X0  ”  100,  M  **  40,  and 

A 

p  ■  -10  dB,  then  rmse0o  ■  0.03°. 

Now  suppose  that  and  are  independent  uniform  variates  whose 
range  Is  (0,L),  i.e.  the  sensors  are  uniformly  distributed  on  the  square 

A 

{0<x<L,  0<y<L}.  Then  cr2  »  L2/12.  Let  us  compare  the  rmse0o  of 

O 

this  random  array  with  that  of  the  square  lattice  array  whose  M-N* 
sensors  are  at  the  points  { (jd,Id) : j ,1-1 , . . . ,N) .  If  the  length  of  the 
square's  sides  Is  L,  then  the  sensor  spacing  is  d  ■  L/(N-1). 

From  (2),  (4),  and  (5),  we  only  have  to  compare  M-*£(xk_x)2  » 
M“1r(yit-y)2  with  o2. 


Since 
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M-1  l  (xk-^)2 

k-1 


M_1d2N  l  (j-j)2 
j-1 


(7) 


,  L2  N+l 

-  d2(N-l)(N+l)/12  -  - 

12  N-l 

«  L2/ 12  -  a2, 

expression  (6)  holds  for  the  square  lattice  array.  The  approximate 
raise  of  the  maximum  likelihood  bearing  estimator  for  a  uniform 

A 

random  array  on  a  square  is  equal  to  the  approximate  rmse60  for  a 
uniformly  spaced  lattice  array  on  the  same  square. 

2.  Three-Dimensional  Random  Arrays 

For  a  given  coordinate  system,  let  xjc  *  (xk»yk»zk^'  denote  the 
vector  location  of  the  kth  sensor  in  a  three-dimensional  array.  Let  0O 
denote  the  azimuth  angle  of  propagation  with  respect  to  the  x  axis,  and 
let  aQ  denote  the  elevation  angle  with  respect  to  the  z  axis.  Thus  the 
signal  at  the  kth  sensor  is 

8(t,Xk)  -  Aexp[i(u>0t  -  <.'xk)]  +  e(t,Xk), 

where  jc '  ■  (xx,<y,<z)  is  the  vector  of  wavenumber  components 
<x  •  (2x/Xo)cos0o,  <y  -  (2irMo)8in0o,  and  <z  -  (2n/X0)coso0. 

The  correspondence  between  beamforming  and  frequency-wavenumber 
processing  holds  in  three  dimensions.  The  ML  estimators  of  the 
wavenumber  components  are  the  kx,  ,  and  <z  that  maximize 


T 
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N  M 

I  l  l 


j-1  k-1 


3^tj»xk*yk»zk)exP^i(£'2Ek  "  “o^)) 


2 


(8) 


where  N  is  the  number  of  simultaneous  discrete-time  observations  of  the  M 

A  A 

channels.-*  The  rms  errors  of  <%  and  <y  are  approximated  by  (2),  and 
rmseKz  a  .  Once  again,  the  ML  estimator  of  the  source 

A  A  A 

bearing  is  0Q  »  tan~^(<y/Kx) ,  and  thus  (5)  holds  for  a  totally  random 
three-dimensional  array  of  M  sensors. 

3.  Estimating  the  Wavenumber  Spectrum 

Consider  the  problem  of  estimating  the  frequency-wavenumber  spectrum 
of  the  ambient,  zero  mean,  Gaussian  noise  field  around  a  random  array. 
Since  an  n-dimensional  array  is  not  much  harder  to  analyze  than  a  linear 
array,  let  x^  *  (x^i , . . . ,x^n) '  denote  the  vector  position  of  the  kth 
sensor  with  respect  to  a  fixed  coordinate  system.  Assume  that  the  x^  are 
realizations  of  independent  random  vectors  {Xj^  *  (X^, . . .  ,X^n) '}  that 
have  a  common  continuous  multivariate  density  f (x) .  Rotate  the 
coordinate  system  so  that  the  covariance  matrix  of  X^  is  diagonal,  and 
for  simplicity  let  <Jj2  «  ...  -  an2  =.  a2^  a2  the  variance  Qf  each 

X^  after  rotation. 

Let  e(t,x)  be  the  noise  at  point  x  at  time  t.  If  the  noise  field  is 
stationary  in  t  and  jc,  the  covariance  function  ce(t,£)  =*  Ee (t+r ,x_+jOe (t ,x) 
is  independent  of  t  and  x.  The  frequency -wavenumber  spectrum  is  defined 
as 


Se(w,jO  -  /c€(T,£)exp[i(K/£  -  U>T)]d£, 


(9) 


1 
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assuming  that  ce  Is  absolutely  integrable.  The  power  spectrum  of  the 
noise  is  Se(u>,0). 

Assuming  that  the  channels  are  sampled  at  times  tj*jA  for 

j“0,...,N-l,  define  the  discrete  Fourier  transform 
N-l 

{eQtjc)  **  e(jA)exp(-iwJA) :  k  »  1,...,M}.  If  Sc(u>,0)  is  bandlimited 

at  u/A,  then  N''^e|c(xjc)  P  a  A“*Se(<»>*0)  for  large  N.^  Let  us  work  with  the 
e(xfc)  to  obtain  an  estimator  of  Se(jc,m)  for  a  given  u>,  which  will  be 
denoted  SeOO  to  simplify  notation.  The  properties  of  the  estimator 
depend  on  the  following  theorem. 


Theorem.  Define  the  n-dimensional  Fourier  transform,^ 

UOO  “  Efe*!  e(xk)exp(i£,xk) .  Assume  that  D(o)  -  / f2(x)dx  -  0(o-n)  and 
whence  *  0,  |$0O|  <  co-11  for  some  constant  c,  where  $(£)  -  Eexp(ijc ’Xj^  is 
the  characteristic  function  of  Xj^.  These  assumptions  hold  for  the 
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lim  (DM2)-1(2n  )_n/  |R(k)  |  2  dK  -  D— ^ < 2n )“ n/  | <J> ( *c )  1 2d< 

M+oo  ~  ~ 

-  D*1Jf2(x)dx  -  1.  (11) 

From  (10),  lim^  (DM2)-1  |R(0)  |2  =  D-1 14>(0)  | 2  -  D-1  =  0(an).  Thus  (11) 
implies  that  as  M  and  o-m»,  (DM2)“1|r(k) |2+S0<),  a  Dirac  delta  function.  If 

£l*£2» 

(DM2)”1R(Ki)R*(52>  -  0(a“n)  +  0p(M"l'2).  (12) 

These  limit  results  are  used  as  follows: 

M  M 

E[U(kj)U*(£2)]  “I  l  ce(*j  “  x^expUOfj'Xj  "  £2'Bc^  (13) 

j-1  k=l 

M  M 

*  (2tt )“n  l  l  /Sc(v)exp[-iv(x.j  -  Xj^)  JexpliUj'jtj  -  ldi 

j-1  k-1 

from  the  inverse  of  (9).  Gathering  terms, 

E[U(kj)U*(<2)]  -  (2it)"ii/r(£1  “  v)R*(<2  “  v>se(v)dv.  (14) 

Thus  from  the  above  limits  and  (14),  limM>a+c»(DM2)_1E|U(j<)  |2  = 

(  2ir  )~nJ 6 (jc-v) Se  (v)d_v  -  Se  (<_) . 

If  *  £2*  then  limM>0>OB(DM2)"^E[U(£i)U*(£2)  1  =  0  from  (12).  Thus 
U(< ! )  and  U(<y)  are  asymptotically  uncorrelated.  For  finite  M<<o2n,  the 
correlation  is  Op(M”l/2). 

This  theorem  provides  a  basis  for  estimating  Se(j<).  One  method  is  to 
divide  the  (time)  sample  into  J  segments  of  successive  observations, 

Nj  -  N/J,  and  compute  U(j<)  for  each  segment.  These  UjOO's  will  be 
approximately  uncorrelated  if  Nj  is  large.  Thus  from  the  theorem, 


9 

Uj (k)  have  a  complex  Gaussian  distribution  for  each  j  (the  noise  is 
Gaussian)  ,  2(DM2)“1  | Uj  (<_)  | ^/Se(jO  is  approximately  chi-squared  with  two 

A 

degrees  of  freedom  and  thus  the  variance  of  Se(jO  is  approximately 

J_1S2(k). 
e  — 

4.  A  Planar  Array  Example 

Continuing  with  the  vector  notation,  suppose  that  the  sensors  are 
uniformly  distributed  on  the  square  {-L/2  <  x^  <  L/2,  -  L/2  <  X2  <  L/2} . 
Thus  f(x)  -  1/L2  for  x  in  the  square,  =.  ,  a2  l2/12,  and 

D  *  / f2(x)dx  *»  L“2.  The  assumptions  for  the  theorem  hold  since  D  ■  0(o“2) 
and  $(0  -  4(x1K2L2)“l8in(KiL/2)sin(K2L/2)  =  0(o-2) .  Thus 
(L/M)2E|UU)  |2  -  S£0O  for  large  M  and  L  in  this  example.  The  estimator 
of  S£0O  is  then  (L/M)2J~1Ej1.1 1  Uj (<) j 2  using  the  time  segmentation  method. 


*Thi8  work  was  supported  by  the  Office  of  Naval  Research  (Statistics  and 
Probability  Program)  under  contract. 
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